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Finite Element Stability Analysis for Coupled 
Rotor and Support Systems 

Pj>.rt III of Final Report Under 
Contract NAS2-7613 

The conventional way of assessing the dynamic stability *f coupled 
rotor airframe systems is to first determine a few undamped purely 
structural blade and airframe modes, to compute the aerodynamic reaction 
for each structural mode, and to consider an aeroelastic mode as a 
superposition of several structural modes, often of only two structural 
modes. The question then arises, how many structural modes are required for 
an adequate representation of an aeroelastic mode, a question that is not 
always easy to answer. The finite element analysis method investigated 
here avoids this difficulty and allows a direct computation of the 
aeroelastic modes. The characteristic equation of the total system is 
obtained by receptance or impedance matching at the interface between rotor 
and airframe. The method is studied for a case of hingeless rotor cyclic 
blade flap-bending motions in hovering coupled to rotor support rolling and 
pitching motions. Of the many aeroelastic rotor modes only two - the 
progressing first and the regressing second flap-bending modes - strongly 
cx/Uple with the rotor support modes and can become unstable for a range of 
rotor support stiffnesses. 



Preface to Final Report under Contract NAS2-7613 


Work under Contract NAS2-7613 started on July 1, 1973. The contract 
was originally awarded for a 3 year period. 

Due to the slower than anticipated progress of the experimental work, 
not all research goals had been achieved by 30 June 1976. Since less than 
the anticipated cost for personnel and equipment had been spent) the 
research contract was extended by a year without Increase In funding. 

The research goals as stated In the contract were: 

(a) Assess analytically the effects of fuselage motions on stability 
and random response. The problem is to develop an adequate but not 
overly complex flight dynamics analytical model and to study the 
effects of structural and electronic feedback, particularly for 
hingeless rotors. 

(b) Study by computer and hardware experiments thi* feasibility of ade- 
quate perturbation models from non-linear trim conditions. The 
problem Is to extract an adequate linear perturbation model for the 
purpose of stability and random motion studies. The extraction Is 
to be performed on the basis of transient responses obtained either 
by computed time histories or by model tests. 

(c) Extend the experimental methods to assess rotor wake-blade 
interactions by using a 4-bladed rotor model with the capability 
of progressing and regressing blade pitch excitation (cyclic pitch 
stirring), by using a 4-bladed rotor model with hub tilt stirring, 
and by testing rotor models in sinusoidal up or side flow. 
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Including the final report, 10 reports under Contract NAS2-7613 
have been submitted. They are listed as P. 1 to P. 10 at the end of 
the Preface. P. 1 and P. 10 pertain to research goal (a). P. 2, P. A, 

P. 6, P. 7, P. 8, P. 9, pertain to research goal (b). P. 3 and P. 5 
pertain to research goal (c). The latter Is not as yet con^lete since 
neither hub tilt stirring nor testing Is slnust)ldal up or side flow 
has been performed. While P. 10 describes only work done during FY 1977, 

P. 8 and P. 9 combine both FY 1977 work results and summaries of earlier 
results, so that the three parts of the Final Report can be read without 
recourse to the earlier reports. P. 8 Includes much new material not 
available when the preceding Yearly Report P. 7 was written. The 
experimental data of P. 9 have all been obtained in FY 77. 

So far 3 publications came out of the research under Contract NAS2-7613. 
They are listed as P. 11, P. 12, P. 13. 

List of Reports and Papers 
under Contract NAS 2- 76 13 

P 1. Hohenemser, K. H. and Yin, S. K., "Methods Studies Toward Simplified 
Rotor-Body Dynamics", Part I of First Yearly Report under Contract 
NAS2-7613, June 1974. 

P 2. Hohenemser, K. H. and Yin, S. K. , "Computer Experiments in Pre- 
paration of System Identification from Transient Rotor Model 
Tests", Part II of First Yearly Report under Contract NAS2-7613, 

June 1974. 

P 3. Hohenemser, K. H. and Crews, S. T. , "Experiments with a Four-Bladed 
Cyclic Pitch Stirring Ifodei Rotor", Part III of First Yearly Report 
under Contract NAS2-7613. 

P 4. Hohenemser, K. H. , Banerjee, D. and Yin, S. K., "Methods Studies 
on System Identification from Transient Rotor Tests", Part I of 
Second Yearly Report under Contract NAS2-7613, June 1975. 
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P 5. Hohenenser, K. H. and Crews, S. T., ’'Additional Experiments with a 
Four-Bladed Cyclic Pitch Stirring Model Rotor". Part II of 
Second Yearly Report under Contract NAS2-7613, June 1975. 

P 6. Hohenemser, K. H. , Banerjaa, D. and Yin, S. K. , "Rotor Dynamic 

State and Parameter Identification from Simulated Forward Flight 
Transients", Fart I of Third Yearly Report under Contract 
NAS2-7613, June 1976. 

P 7. Hohenemser, K. H. and Crews, S. T., "Rotor Dynamic State and 

Parameter Identification from Hovering Transients", Part II of 
Third Yearly Report under Contract NAS2-7613, June 1976. 

P 8. Hohenemser, K. H. and Crews, S. T. , "Unsteady Hovering Rotor 

Wake Parameters Identified from Dynamic Model Tests", Part I of 
Final Report under Contract NAS2-7613, June 1977. 

P 9. Hohenemser, K. H. and Banerjee, D. , "Application of System 

Identification to Analytic Rotor Modeling from Simulated and 
Wind Tunnel Dynamic Test Data", Part II of Final Report under 
Contract NAS2-7613, June 1977. 
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Nomenclature 

blade flap-bending stiffness 

EI/EIq flap-bending stiffness referred to root value 

tremsfer matrix across massless blade element 

transfer matrix across point mass 

rotor support moment of inertia about rotor center 

blade moment of inertia about rotor center 

non-dimensional blade flap-bending moment (unit PqR^Q^) 

number of finite elements per blade 

polynomials of X 

rotor radius 

non-dimensional blade shear force (unit 

non-dimensional blade centrifugal tension force (unit p^R^fl^) 

state vector 

blade airfoil lift slope 

number of blades in rotor 

non-dimensional blade chord (unit R) 

aerodynamic damping coefficient for ith blade element 

aerodynamic non-dimensional force at blade station i 
(unit p^R^n^) 

aerodynamic coefficient for blade pitch angle 
non-dimensional length of blade element (unit R) 
non-dimensional point mass at blade station i (unit PqR) 

(T/q eT)^^^ blade tension force parameter 
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Nomenclature 

q 

t 

Up. Ut 

X 

y 

z 

0 

L>] 

'<'k 

n 

a 

1 

P 

Pc 

<u 

Subscripts 

I. II 
i 

12 . 13 ... 


( cent ' ) 

Ei^/PoR'^n^ non-dimensional blade root flap-bending stiffness 

non-dimensional distance of blade station i from rotor 
center (unit R) 

non-dimensional time (unit 1/R) 

non-dimensional normal and tangential velocity at blade 
element (unit nR) 

non-dimensional distance from outer end of blade 
element inboard (unit R) 

non-dimensional blade deflection, positive up (uni R) 
a single blade variable 

blade pitch angle, positive nose-up 

transfer matrix 

azimuth angle of kth blade 

rotor angular speed 

hub tilt angle 

real part of coupled frequency (unit ft) 

E i iw. non-dimensional complex valued frequency, (unit n) 
non-dimensional air density (unit Pq/R^) 
blade mass per unit length at blade root 
circular frequency (unit J2) 


multiblade coordinates, forward and left respectively 
value at ith blade element, beginning at blade tip 
denotes elements of a determinant or matrix 
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Nomenclature (cont*) 

1; ** Nfl denotes values at blade tip and blade root 

I, R Imagir.ary and real parts 

M, a, 0 polynomials multiplied by «» © 

1, r left, right of mass m^ 

time differentiation 
length differentiation 
truncated polynomials 
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Introduction 

For the dynamic design of lifting or propeller rotors one usually 
begins by studying potential single blade dynamic insted>ilities. They 
can be caused by the effects of coupling between blade flap-bending, 
lag-bending, and torsion^ by the periodicity of the aerodynamic stiffness 
and damping coefficients, or by unsteady aerodynamics phenoiAena. Even 
for stable single blade dynamics, potential instabilities can still 
exist, caused by the coupling of tnp blades with each other, with the 
rotor support or airframe, and with the rotor controls, A widely used 
approach to the problem of rotor dynamic stability, loads and vibrations 
is to first compute a few purely structural blade normal modes, using 
for example the methods of references 1 and 2, and then to consider the 
aeroelastic modes as truncated series of the structural normal modes , 

The truncation error is sensitive to the details of handling the 
aerodynamic loads, see for example references 3 and 4. The truncated 
normal structural mode expansion is either used as the basis of non-linear 
"global” computer programs as for example Rexor or MOSTAB., or it is used in 
linear eigenvalue programs as for example in reference 5, where two 
flapwise bending modes, two chordwise bending modes, one elastic torsion 
mode, and a few airframe and control modes are considered. 



The problem with the tindamped structural normal mode analysis is 
that in particular the low flap-bending aeroelasxic modes are 
substantially affected by aerodynamics so that sizeable errors can be 
expected when only a few purely structural normal flap-bending modes 
are used. In a finite element stability anedysis these errors are 
avoided and the aeroelastic modes are obtained directly. A finite 
element eigenvalue analysis results in a high order characteristic 
equation and requires a high precision computer program. The concept 
to be studied here is to perform a separate finite element eigenvalue 
analysis for the rotor and for the airframe, and then to couple both 
systems with the help of kinematic and equalibrium equations tedcen at 
the interface, A simple rotor and support representation is selected to 
gain some initial experience with the finite element stability analysis. 

Modeling of Rotor and Support System 
A schematic view of the rotor and support system used for the 
analysis is shown in Fig. 1. The blades are assumed rigid in chordwise 
bending and torsion, but flexible in flap-bending. They are rigidly 
connected to the hub so that the slope of the deflection, line relative 
to the hub is constant and equal to the built-in coning angle. Since 
blade deflections and coning angle are small and only out-of-plane blade 
bending is considered, the built-in coning angle has no effect on rotor 
dynamics. The rotor shaft is assumed to be rigid and connected to a 
rigid housing that is supported by a focussing mount with focus on the 
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rotor center. Thus the hub is rigidly restrained against vertical and 
sidewise motions. The elastic restraints in pitch and roll of the hub 
are indicated in Fig. 1 by the horizontal springs attached to the shaft 
housing. There also is a gravitational restraint (pendulum effect). The 
restraining springs may have different stiffness in pitch or roll. 

Fig. 2 shows the hub witn a blade cross section in a position of 
angular deflection in pitch. It is assumed that the blades perfon., a 
uniform rotation about the vertical axis Z-Z. In actuality the 
rotational motion is non-uniform in a complex way depending on 
chordwis^ blade-bending and drive system dynamics. The vertical blade 
deflections are measured from the horizontal reference plane through 
the rotor center, indicated in Fig. 2 by X-X. Aerodynamic forces on the 
blades are produced by their vertical motions v;ith respect to the 
reference v -=me (aerodynamic damping) and also by the blade pitch angle 
ch6i*.ges ‘..■'th respect to the reference plane, as seen in Fig. 2. It is 
assumed that the vertical blade motions and pitch angle changes do not 
produce a cha..ge in vertical inflow through the reference plane, an 
assumption approximately satisfied for the higher frequency regime that 
will be here of primary interest. Blade pitch controls are assumed to 
be rigid so that the blade pitch angle is equal to the hub tilting angle 
about the blade axis, indicated in Fig. 2 by Cj = ©j. 

Blaaes with constant chord and uniform mass and flap-bending stiffness 
are assumed, though the method of analysis is suitable also for arbitrary 
blade planform and for non-uniform mass and stiffness distribution. The 
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uniform blade flap-bending stiffness extends to the rotor center. In 
actual hingeless rotor designs the hub is very stiff and the root 
section of the blade is no re flexible then the rest of the blade. 

Extending the uniform blade flap-bending flexibility to the rotor center 
results, however, in first and second bending mode natural frequencies 
typical of actual hingeless rotors. 

The analysis is performed in a non-dimensional way with length 
unit R, mass uiikit and time xinit X/U, The force unit then becomes 
PoR^n^i the unit of angular inertia becomes PqR^* and the unit of 
bending stif^ess becomes PoR'^< 2^* Forgiven mass and flap-bending 
distribution the rotor blades are described by the non-dimensional blade 
root bending stiffness q and by the blade Lock number y. The rotor 
support is described by the ratio of blade moment of inertia over support 
moment of inertia, both taken about the rotor center, Ib/I» and further 
by the two non-dimensional support frequencies without blades Uj and (Ojj. 
Thus the rotor and support system is uniquely defined by 5 parameters: 
q, Y, Ijj/I, Wj, ujj. The parameter q"^/^ is proportional to the rotor 
angular speed 0 and can be used as rotor speed parameter. In the numerical 
examples q“*/2- £g selected to simulate the dynamics of actual 
hingeless rotors. The remaining parameters are then varied within wide 
limits . 


Method of Analysis 

Reference 1 describes an analysis for determining the undamped 
structural single blade bending modes and frequencies with the help of 
blade element transfer matrices. The natural frequencies are obtained 
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by establishing through trial and error the zeros of the characteristic 
determinant, a method that is not feasible for un conservative systems. 

The following four steps are taken in extension of reference 1. 

1. Aerodynamic terms are included in the blade transfer matrix. 

2. Relations between blade root state variables are derived Including 
a dynamld blade pitch feedback term. 

3. Single blade polynomials are transformed Into multiblade polynomials. 

4. The total system characteristic equation is derived by receptance or 
Impedance matching. 

In a later section an intermediate step will be inserted between Steps 
2 and 3 to reduce the order of the characteristic equation. 

Blade Element Transfer Matrix 

Fig. 3 shows a blade element between station i (outer end of element) 
and station i + 1 (inner end of element). The centrifugal force is 
assumed to act in the undeformed position of the blade and produces a 
bending moment T^y^. The moment is thus fictitious and + Tj^y^ 
is the actual bending moment acting on the outer end of the blade element. 
All quantities are non-dimensional with the units listed before. The 
inertia force -yim^ = -X^m^y^ and the aerodynamic damping force 
-y C£ = -X C£y£ are added to the shear force at the point mass m-^. 

Also added is the aerodynamic force from dynamic blade pitch 0, expressed 
as g^ 0. Since the blade is rigid in torsion, 0 is constant along the 
radius. Dynamic pitch 0 is caused by hub motions, see Fig. 2. Denoting 
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by 1 and r quantities to the left and right of the point mass m^^, 
we have 

Sii « Sij, - (x2 + X CiYi) + g£ 0 (1) 

The change of the fictitious moment across the point mass is 

since this incremental centrifugal force is assumed to be located in the 

undeformed blade position rather than at the mass m^^. Thus 

^il * ”ir ■ 

Over the length of the blade element 1^ bending stiffness Efq and 
centrifugal force T are constant. The deflection y is determined 
from the differential equation of beam bending 

y" Elq = Sjf + (M t Ty) (3) 

where x is measured from the outer end of the blade element inboard. 

S and M are shear force and fictitious bending moment to the left of a 
point mass m^. The slope y'(x) is positive when opposite to y*£ 
shown in Fig. 1. The solution to Eq. (3) is 

y(x) = y(0)cosh (px) + (y' (0)/p)sirih (px) (4) 

+ (M/T)[cosh (px) - 1] + (S/T)[(l/p)sirih (px) - x] 

The derivative is: 

y'(x) s y(0) p sinh (px) + y'(0) cosh (px) 


+ (M/T) p sinh (px) + S/T[cosh (px) - 1] 


(5) 
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where p = (T/q 


( 6 ) 


The fictitious bending moment is 

M(x) = M + Sx (7) 

The shear force is constant 

S(x) s S (8) 

Denoting the state vector by 

(9) 

one can now write the change of the state vector along the blade element 
as a product of two matrices 



^i+1 ^ ^i 


( 10 ) 


The matrix [F] gives the change in state vector from the right side of 
the mass m^ to the left side and expresses Eqs. (1) and (2) together 
with the continuity relations * ^1 ” ^r* ®1 ~ 


1 0 
0 1 
0 0 
0 0 
0 0 


0 

0 

1 

0 

0 


-d^mi + Xc^) 


-m.ri 


0 

1 

0 



0 

0 

0 

1 


CF] = 


( 11 ) 
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The Mtrlx [E] gives the chsnge in state vector from the left side of the 
mass mi to the right side of the mass m^^^. It expresses E,s. (4) to 

(8) considering that y| = -y'(o). y'^^s = Oj. The ,«trix 

[E] then is 



where 

^31 " "^42 * ”(l/Ti)Ccosh (p.lj) - 1] 
E32 = -d/T.) ainh (p£lj^) 

^33 * ^44 ‘ 

E34 = p. sinh (P£l£> 

E41 = (l/T^)(l/p^) sinh (p^l^) 

E43 = -d/p^) sinh (P£l£) 

Pi = ^^i/q 


( 12 ) 


(13) 


In comparing these expressions with those given in reference 1 one should 
note that A is defined differently leading to the opposite sign of A^ 
in Eq. ( 11 ). Furthermore the unit of circular frequency in reference 1 Is 
"hile it is here 0. 
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Denoting the matrix product [E][F] = [(j)], Eq. (10) can be written 
in the form 




(14) 


For N masses one obtains by successive transfer matrix multiplication, 
beginning at the blade tip, for the root quantities the following relation 


— . .. 







— — 


— - 

^N+1 


♦ll 

<^12 

*13 

'^14 




♦l5 

«Ntl 


^21 

*22 

*23 

*^24 




*^25 

yi+i 


*31 

*32 

*33 

*^34 




^35 

^N+l 


^^41 

* 1^2 

^43 

^4 




^45 


. 









(15) 


All elements of are polynomials of X. Eq, (15) completes step one. 
Before proceeding further, the aerodynamic coefficients c^. and must 
be determined. 


Aerodynamic Coefficients 

As before, x is mea ored from the outer end of the blade element 
toward the rotor center. Non-dimensional tangential and normal flow 
velocities at the blade element u^ and Up are, if dynamic inflow is 
ignored, and if Up is averaged over 1^: 

• Up = - (i/2)(yp + 


Up = rj - X 


(16) 
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The non-dimensional aerodynamic force acting on the blade element of 
length 1^ is: 


f, - (1/2) 



p ac(0 u 


2 

T 


+ Up Uj.)dx 


o 

Introducing the Lock number for a uniform blade which takes in our 
non-dimensional units the form y = 3 P and inserting Eq. (16) 
into Eq. (17): 


(1/6 )y 



{0 (r^-x)^ 


(l/2)(y^+y^^j^)(r^-x)}dx 


( 17 ) 


(18) 


Performing the integrations: 

f^ = (1/6 )y 0(r? 1^ - r^ 1? + l?/3) 

- (1/12)y (y. + >i+iHr.l. - l?/ 2 ) 


(19) 


One half of the first term is assumed to act at statioai, the other half 
at station i+1. Of the second term the force with factor y. is assumed to 
act at station i, the force with factor at station i+1. With this 

assumed distribution of the aerodynamic forces the total aerodynamic 
force at station i is: 


^Ai = ^1-1 - -^1-1 ^i-1 


. li - if f if /3)-y,(r,_^ li.l - ll,/2 


" *-1 If - 


( 20 ) 


The factor of 


0 is equal to , the factor of y^ is equal to c.. 



14 


Blade Root Relations 

Equation (15) relates the variables at the blade tip to those at the 
blade root. We need, however, a relation between blade root variables. 
Two variables at the tip are known, since shear force and bending 
moment are zero. The centrifugal force is also zero at the tip, thus 
the fictitious bending moment is zero. There is no vertical motion at 
the rotor center, so that 


The two remaining tip variables y^, will be expressed in blade root 
variables, as follows 

Inserting Eq. (21) into Eq. (15) we have 




r* 

— 




r “7 

^N+1 


*^13 



» 

yi 

t 

♦u 

^N+1 


CO 

CM 

-o- 

CM 

-e- 


^1 


■^25 

mm mm 


L. 

-J 


— - 


— - 


%+l 


*33 

♦34 


f 

yi 

+ 

1 m 
CO 

I ^ - - 

0 


*43 

♦ 44 


^1 


"^45 



— 



— 1 


— - 


0 


0 


( 22 ) 


( 23 ) 


From Eq. 



( 24 ) 
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We now substitute Eq. (24) into Eq, (22); 


- 


1- 

-1 


- 

- 

f-1 

f 


- -| 



r “T 

^N+1 

s 

^13 

^14 


^33 

*34 

( 



*35 


+ 

in 

”n+i 


♦23 

♦24 


♦43 

*44 

\ 

0 


*45 

/ 


♦25J 



Im 

«J 







— - 



L. 


( 25 ) 


After performing the required manipulations, one obtains for the second 
equation 


yw+i ® 

♦23 *33 ^43 

♦24 ♦sq ^44 

^25 ^35 ^45 

(26) 

Equation (26) can only be used for a single blade characteristic 
equation, if the dynamic pitch angle 0 is either zero or related to a 
root variable. For example a hinged blade without pitch-flap coupling, 

* 0 s 0, gives as characteristic equation 

Pa(A) = 0 (27) 

A hinged blade with a pitch-flap coupling defined by 0 = ^y^+l* S^ves 
the characteristic equation 


where * 


*33 

*34 


*23 

*24 



. Pa(X) = 



*43 

*44 


CO 

-G- 

*44 


, P0(X) = 


Pa(X) + k Pg(X) = 0 


(28) 
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Finally a cantilever blade, where = 0*0 gives the characteristic 

equation 

Pj^(X) = 0 

This completes step 2 of the analysis. 

Multi Blade Transformations 

The relations between single blade and multi blade coordinates are, 
when only cyclic terms are retained, see reference 6, 

‘yN+i>k ' “i ♦k ♦ “ii *k 

<"nn>k ' “i “= h * "ii *k 

Oj, Mj, 0j refer respectively to nose down rotor tilting angle, nose 
down rotor moment on its support, and nose down cyclic pitch angle. 

Ojj, Mj, 0JJ refer to left rotor tilting angle, left rotor moment on 
its support, and left cyclic pitch angle. Equations (30) are easily 
inverted. 

For example 

b 

S '>N+l>k ♦k 

k=l 


Oj = (2/b) 


(31) 
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with b the number of blades of the rotor. Eqs. (30) and (31) are 
valid for b ^ 3. Equation (26) is now transformed into relations 
between multiblade coordinates. The question is how to transform a 
polynomial of X for a single blade coordinate into polynomials of 
X for multiblade coordinates. Each factor of X corresponds to a 
differentiation. Assume an arbitrary single blade variable z and its 
multiblade counterpart Zj, z^j related by 

z - Zj cos t + Zjj sin t 

whereby the azimuth angle has been written in the form of the non- 
dimensional time t in which the time of one rotor revolution is 
equal to 2ir. Differentiating Eq. (32) once: 

z = (Zj + z^^)cos t + (~Zj + Zjj)sin t 
Differentiating a second timet 

z = Czj - Zj + 2 Zjj) cost - ( -2 z£ + Zjj - Zjj)sin t 

Replacing the differential quotients by factors of X E<p. (33) and ( 
Xz = (XZj + Zjj)cos t + (-Zj + Xzjj)sin t 

X^z - (Zj(X^ - 1) + 2X Zjj)cos t + (’2\ Zj + Zjj(X^-l) )sin t 


(32) 

(33) 

(34) 

34) read 

(35) 

(36) 


Let us replace in the polynomial on the left hand side X by (X+J) and let 
us then form a matrix that has as diagonals the real parts of the left 
hand side polynomial and as off diagonal terms the plus and minus 
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imaginary parts of the polynomial, plus being used for the first row 
For Eq. (35) this matrix reads 

X 1 

-1 X 

For Eq. (36) this matrix reads, since (X + i)^ = - 1 + 21 X 


X2 - 1 2X 

-2X X2 - 1 

Post multiplying these matrices by the column LZr/Zrrl one obtains 

1 XX 



(37) 


X2 - 1 2X 
-2 X X2 - 1 




(X2 - 1) 


4 - 


2X 






(38) 


Equations (37) and (38) yield the factors of cos t and sin t in Eqs. (35) 
and (36). In general, if a single blade coordinate z is multiplied by 
^ polynomi il P(X) as in Eq. (26), one obtains the multiblade expressions 
by splitting the polynomial P(X + i) in real and imaginary parts: 

P(X + i) = P]^(X) + i Pj(X) and by writing 
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( 39 ) 


To apply this rule, we write Eq. (26) as 

”n+ 1 * ^N+l + 0 Pe(X+i) 


(40) 


In multiblale coordinates Eq, (26) then becomes, considering Eq, (30) 





(X) 

(X) 


Pai<X) 


P„r(X) 





(41) 


From Fig, 2 there is a cyclic pitch feedback represented by the relation 

®I “ (42) 

Inserting Eq, (42) into Eq, (41) one then obtains the relations between 
the rotor moments and the hub angular 'Reflections, from which either 
rotor receptance or rotor impedance can be conputed. The characteristic 
equation for the rotor alone is obtained by setting = M„ = 0. 

Including the feedback from Eq, (42) the characteristic equation reads 
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PoR 

Pal 


-P0I 

P0R 



+ 



■Pal 

PoR 


■P0R 

1 

o 


This completes step 3 of the analysis. 


( 43 ) 


Impedance or Receptance Matching 

It now only remains to match the rotor support to the rotor. The 
moments transferred from the rotor to the support are 

(b/2) Po , (b/2) Po Mjj (44) 

whereby the number of blades per rotor, b, must be at least 3, see 
reference 6. The factor p^R is required to obtain dimensional 
moments. The time unit is still 1/n, otherwise the expressions (44) 
would also have to be multiplied by n’. For uniform blades 

3 

Po R = 3 Ib 

and the rotor support dynamic equations without support damping are 


0 


““ —1 
“l 






= (3 b/2)(I^^/I) 


0 (*)__^ - x2 


a^_ 


M 

II 


II 


■^II 

— — 


— — 




We can now either insert Oj., from Eq. (46) into Eq. (41) (receptance 

matching) to obtain with Eq. (42) a set of homogeneous equations for Mj, 
Mjj. The coefficient determinanc then represents the characteristic 
polynomial. We can also insert from Eq. (46) into Eq. (41) 
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(impedance matching) and then obtain with Eq. (42) a set of homogeneous 
equations for 


Q 

Q 


11 


21 

^ 9.2 



The total system characteristic eqviation is then given by 


(47) 


^12 


^21 ^22 


= 0 


(48) 


The mode shapes can be easily computed by inserting into Eq. (47) one of 
the eigenvalues obtained from Eq. (48). One then finds that all modes 
are either regressing or progressing in the non-rotating reference system. 
This completes step 4 of the analysis. 

If each blade is represented by N point masses, the characteristic 
equation (48) will be of 4N + 4 order. For an airframe the receptance 
or impedance matrix at the rotor-airframe interface will be less 
simple than Eq. (46). If the airframe is represented by L masses or 
moments of inertia, the order of the system characteristic equation 
will be 4N + 2L. A method to reduce this order without appreciable loss 
in accuracy of the lower eigenvalues will be discussed later. First 
some numerical results of solving for the roots of Eq. (48) will be 
presented. 
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Numerical Examples 

The computations were performed on the IBM* 360/65 computer using 
double precision (16 digits). Single blade computations were made 
for 5 1 8, 10, 15 and 20 elements per blade. The number of 20 elements 
was found to be too high for the 16 digit precision used. Kith 10 
blade elements no computational difficulties were encountered for the 
multiblade analysis, provided that the evaluation of the 3 by 3 
determinant in Eq. (26) was numerically optimized by writing 


l*’*l * ■ '^33^*^24 *45 ■ *44 *25^ *34^*23 *45 " *43 *25^ 


(49) 


*35^*23 *44 ■ *43 *24^ 

whereby the set ♦331 '^34» ^35 represents the column with the highest 
values of the elements. A convenient check for adequate computer 
precision consists of looking at the coefficients of 
^nmax+l^ etc., in Pq, whereby n^iax highest power that should 

occur theoretically in this polynomial. For a single blade with N 

point masses n^^ax * ^ n-2 for P0. If the computer precision is 

, ^ ^max+l ,ntnax+2 

adequate, the coefficients of X , X are several orders 

smaller than those of The accuracy of the computation 

depended on how this determinant was evaluated, see Eq. (49). 

Once the single blade problem could be solved without difficulties, 

the complete coupled system solution posed no further obstacles. The 

computer accuracy depends on the highest eigenvalue considered and does 

not suffer vdien the number of eigenvalues is approximately doubled as for 
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the coupled systenii as long as the highest eigenvalue remains approxi- 
mately the same. Though 10 point masses per blade would have been 
satisfactory, the parametric studies for which selected results will 
be presented here were made with 8 masses per blade. The errors In 
the first 3 blade eigenvalues shown here were found to be lass than one 
to two percent. For 5 blade elements more substantial eigenvalue errors 
occur. For 8 masses per blade the equivalent CPU time to obtain a 
complete set of eigenvalues for the coupled system was 4 seconds, for 
55 cases computed In one run the CPU time was 34 seconds. 

Sl.ngle Blade 

The two parameters that determine the eigenvalues for the uniform 
single blade are the non-dimensional bending stiffness q and the Lock 
number y. As mentioned before, the rotor speed parameter q“^/^ = 18 is 
selected to match the flap-bending frequencies of actual hingeless rotors. 
For the blade Lock number we select the values y - 5 and y = 8 that 
also cover the range of actual hingeless rotors. Though computed for a 
single blade with Eqs. (27) and (29), we give in Table 1 the eigenvalues 
in multiblade form to facilitate the comparison with the coupled system 
eigenvalues to be determined later. The letters R and P refer to 
regressing and progressing modes respectively. The numbers 1, 2, 3 
refer to first, second, third blade flap-bending mode. £ and w are real 
and imaginary part respectively of an eigenvalue. 
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Table 1 

Uncoupled Blade Eigenvalues in Non-Rotating 
Reference System, 


Flap-Bending ' 
Mode 

Hinged Blade 

Cantilever Blade 


Y = 

5 

Y 

= 8 

Y 

= 5 

Y 

= 8 


5 

Ui 

K 

u 

1 5 

u 

K 

U) 

IR 

-.311 

.05 

m 

.13 

-320 

.01 

-.514 

.07 

IP 


1.95 

IBI 

1.87 


2.01 


1.93 

2R 

-.257 

1.57 

1 

-.409 

1.55 

-.265 

1.74 

-.424 

1.72 

2P 


3.57 


3.55 


3.74 


3.72 

3R 

-.229 

3.80 

-.366 

3.79 

-.237 

4.17 

-.378 

4.16 

3P 


5.80 


5.79 


6.17 


6.16 


It is seen that there is little difference between the eigenvalues 
for tb«* hinged and cantilever blade. This is due to the relatively low 
flap-bending stiffness of the blade as expressed in = 18. The 

first natural frequency without aerodynamic damping is 1 for the 
hinged blade and 1.06 for the cantilever blade. The aerodynamic 
damping has the effect of lowering the natural frequency. The effect 
of increasing the Lock number from y = 5 to y = 8 is merely to increase 
the damping of each mode in the ratio of 8/5 and to lower the natural 
frequencies somewhat. The numerical examples for the coupled system will 
be limited to y = 5. When determining the mode shape one finds 
substantial phase shifts along the blade. For example with Y = 5 the 
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first mode of a hinged blade shows between root slope and tip slope a 
phase shift of 12 degrees. For y = 8 the phase shift is still higher. 
Thus there is a substantial difference between the structural mode 
shape and the aeroelastic mode shape. 

Rotor Alone (Inter Blade Coupling) 

Coupling the blades to each other and assuming zero moments at the 
rotor center, the evaluation of Eq. (43) for ^ j_g g^^d y ~ 5 

yields the following eigenvalv.es, whereby the character of a mode as 
regressing (R) or progressing (P) can be determined firom Eq. (47). 

Table 2 

Rotor Alone Eigenvalues for q"*/^ = 18, y ~ 5 


Flap- Bending 
Mode 


l 


u 


IR 

IP 

2R 

2P 

3R 

3P 


0 

-.616 

-.327 

-.193 

-.240 

-.214 


0 

1.96 

1.57 


3.62 


3.87 

5.87 


The values of Table 2 should be compared to the values in Table 1 
for the hinged blade and y - 5, The frequencies are almost the same, 
and the sum of regressing and progressing mode damping are also the same. 
However now the regressing and progressing modes have different damping. 
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With respect to the first regressing mode the frequency and damping are 
now zero, the damping of the progressing mode is doubled. With respect 
to the second and third mode the progressing mode damping is reduced, 
the regressing mode damping is increased. 

Unrestrained Rigid Body 

When omitting the restraining springs shown in Fig. 1 and the 
gravitational effect, so that the rotor is coupled to an inertia, one 
obtains an approximation to the short period pitch and roll modes in 
helicopter hovering flight mechanics. In Eq. (46) we have Wj. = = 0. 

V/e present here a case of a rotor with 4 blades (b = 4), and an 
inertia ratio Ib/I = .2. Eq. (46) could of course easily be written 
for different body inertia in pitch and roll. For _ ^8, y = 5 

one obtains : 

Table 3 

Coupled Rotor-Rigid Body Eigenvalues 


q"l/2 = 18^ Y = 5, b = 4, 1^/1 = .2, = “n = ° 

MODE , i (j 


Predominant 

R 

r- 

-.157 1 

.16 

Body 

P 

f 

-.167 ! 

.16 


IR 

0 

I 

0 


IP 

- . 316 ; 

2.01 

Pr>pHnmi n pn 


1 



2R 


1.74 

Flap-Bending 

2P 

-.265 1 

1 

3.74 

1 

3R 1 

1 

-.237 ; 

4.17 

! 

3P 


6.17 
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The coupling with a rigid body has added two modes and has substantially 
changed the rotor alone inodes of Table 2. The progressing rotor mode 
with first flap-bending has only one half the damping. There are now 
two predominant body modes , one regressing, the other progressing, 
both with about the same frequency. The second and third flap-bending 
modes have not been much changed by the coupling with the rigid body. 

Rotor and Support Coupling 

In evaluating Eq. (46) for a case of a rotor coupled to its support, 

we assume again 4 blades (b * 4), and now stipulate 1^/1 = 5. The 

support frequencies without rotor will be varied from .4 to 5.5. 

The frequencies In pitch and roll will first be assumed equal, followed 
by a case of unequal support frequencies. Table 4 shows the coupled 
rotor and support eigenvalues for a number of support frequencies 
u>I s Figure 4 shows the imaginary parts of the eigenvalues, u. 

Fig. 5 shows their real parts, C. At low support frequencies up to 
«I = Wjj < .8 one can clearly identify a mode as predominantly support 
or rotor mode. The same is true for high values of Wj = > 2.5. 

In between there is strong coupling between support and rotor, and no 

predominance can be established. There is an essentially aperiodic 
mode with very low frequency. It's almost aperiodic decay is quite low 
for small but increases to -.30 at high oij. The first progressing 

rotor flap-bending mode with an uncoupled frequency 1.96 (Table 2) and a 
high danping increases its frequency and loses its damping with stiffer 
support. This mode becomes unstable at = Wjj = 1.5 and only 


Table 4 


Coupled Rotor and S^port Ei^nvalues 
q-1^ = 18. Y = 5. b r 4. Ijj/I = 5 


1 “I = “II 


.4 


.8 

1.2 

1 

1.6 

MODE 

C 


e 

la 


(I) 

C 


Support Rotor 2 


mm 




-.188 

1.44 

-.235 

1.58 

S\J 5 >port Rotor 1 






-.339 

1.60 

-.403 

1.77 

Rotor 1 

R 

-.037 

.001 

-.109 

0 

-.172 

„0C2 

-.216 

.004 

Rotor 1 -► Support 


-.176 

2.02 

-.145 

2.03 

-.071 

2.06 

+.019 

2.19 

Rotor 2 -► Support 


-.206 

1. 85 

-.188 

1.87 

-.141 

1.94 

-.077 

2.10 

Rotor 2 


-.269 

3.76 

-.269 

3.76 

-.269 

3.76 

-.269 

3.76 

Rotor 3 


-.232 

4.21 

-.232 

4.21 

-.232 

4.21 

-.231 

4.21 

Rotor 3 


-.237 

6.19 

-.237 

6.19 

-.237 

6.19 

-.236 

6.19 


a>I = (jjjj 

2.5 


3.5 


4.5 

5.5 

MODE 

5 

u 

5 

b> 


b> 


u 

Support -► Rotor 2 


-.268 

1.69 

-.269 

1.72 

-.268 

1.73 

-.268 

1.73 

Support Rotor 1 

p 

-.398 

1.96 

-.357 

2.00 

-.341 

2.01 

-.333 

2.01 

Rotor 1 




-.290 

.011 

-.301 

.012 

-.307 

.013 

Rotor 1 Support 




-.053 

3.60 

-.012 

4.65 

-.024 

5.56 

Rotor 2 -► Support 


-.028 

2.76 

-.049 

3.60 

-.037 

4.74 

-.013 

5.63 

Rotor 2 


-.267 

3.77 

-.199 

3.83 

-.251 

3.72 

-.261 

3.73 

Rotor 3 




-.197 

4.28 

-.206 

4.07 

-.232 

4.13 

Rotor 3 


HBOS 


-.236 

6.19 

-.233 

6.20 

-.216 

6.24 
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stabilizes again at **>j ~ become a predominant 

support mcde. The second regressing rotor flap bending mode with an 
uncoupled frequency of 1.57 (Table 2) and a high damping follows 
essentially the same trend, however without becoming actually unstable. 
The progressing second flap~bendlng mode with an uncoupled frequency 
of 3.6 (Tjable 2) emd a damping of -.19 Increases frequency and 
damping due to coupling with the support. The third flap-bending modes 
are not much affected by the coupling and correspond to the cantilever 
cases of Table 1. The change in type of mode when increasing support 
stiffness is indicated in Table 4 by an arrow. For example 
Support -*■ Rotor 2 means that for low we have predominantly a 

support mode, at high Uj we have predominantly a second flap- 

bending rotor mode, both regressing. 

The question has often been raised whether airframe modes can be 
damped by coupling with the rotor. Fig. 5 shows for our case that 
the support damping from the rotor is very small for support frequencies 
above 4.5, and that one must avoid support frequencies between 1.5 
and 3.0 where a coupled mode is unstable. However for low support 
natural frequencies below 1.2 the rotor provides very good damping to 
the support. Obviously the result of Fig. 5 cannot be generalized, 
since it pertains only to the selected combinations of rotor and support 
parameters . 

In some studies, it was found that unequal stiffness of the rotor 
support in pitch and roll can alleviate instabilities. In the case of 
0 )^ - 1.6, Ujj - 4.8, the following eigenvalues were obtained. 
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Table 5 

Coupled Rotor and Support Eigenvalues 
= 18, Y * 5, b = 4, Ij^/I = 5, Wj=1.6, Wjj=4.6 


Mode 

1 

« 1 

fa) 

1 

Support R 

-.266 i 

1 

1 

1.63 

i 

P 1 

-.018 

4,96 

R 

-.216,-, 303 

1 

0 

P 

i 

-.357 

1.90 

R 

-.027 

2.15 

Rotor 2 

P 

-.262 

! 3.74 

1 

R 

-.226 

1 

4.16 

Rotor 3 

P 

1 

-.234 

i 

I 6.20 


As compared to the .of = 1.6 in Table 4 the instability is 

removed, though the second regressing rotor mode has now almost no 
damping. 

Reducing the Order of the Characteristic Equation 

Although for the simple case assumed her there was no difficulty 
in solving the characteristic equation, the introduction of more 
degrees of freedom for the blades (in-plane bending and torsion) and 
of a much higher order for the airframe receptance or impedance 
polynomials may lead to difficulties. It has been suggested in 
reference 7 to reduce the order of the characteristic equation of 
coupled large systems by modalizing and truncating the component 
receptance matrices. It is, therefore, of interest to find out 
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whether the method suggested in reference 7 can result in a lower order 
characteristic equation without an appreciable reduction in accuracy for 
the lower eigenvalues. Let xis write Eq. (26) in the form 

”n+1 ^ ® P0<3^)/Po<^> <50) 

Wo now determine the roots of 

Pa(X) » 0 (51) 

which according to Eq. (27) are the eigenvalues for the single blade that 
is hinged at the root and experiences no feedback from pitch angle 0 . 

We then perform a partial fraction expansion of the two polynomials 
Pj^(X)/Pa(X) and P 0 (X)/Pa(X) with respect to the roots of Eq. (51). 

Since the roots occur in conjugate complex pairs, we have for example 

P„(X)/P,(X) = Co . . •••• (S2) 

and a similar expression for Pg(X)/Pa(X) . This partial fraction 
expansion is now truncated, and Eq. (26) is rewritten as 

”n+ 1 " ^N+l + 0 P 0 (^) (53) 

where the new polynomials ^j^(X), Fjj(X), ^©(X) are of lower order than 
the original polynomials. From here on the analysis follows the same 
procedure as described before leading to the total system characteristic 
equation (40) which now is of lower order. 

Using the case = Ujj = 1.2 in Table 4 and omitting in Eq. (52) the 
last 8 terms, one obtains the following total system eigenvalues. 
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Table 5 

Coupled Potor and Support Eigenvalues 
with Truncated Polynomials 

- 19^ Y = 5, b = 4, Ij^/1 = 5, = Wjj = 1.2 


MODE 

i 

5 

u 

Support ■> Rotor 2 

1 

R 

-.200 

1.44 

Support -► Rotor 1 

P 

-.384 

1.62 

Rotor 1 

R 

-.158 

.008 

Rotor 1 -► Support 

P 

-.038 

2.11 

Rotor 2 -► Support 

R 

-.139 

2.04 

Rotor 2 

P 

-.271 

3.80 


R 

-.232 

4. 34 

Rotor 3 

P 

-.239 

6.31 


In comparing these numbers with those in Table 4 for oj^ = = 1.2 

it is seen that the truncation method has resulted in pood approximations. 
It can be shown that these approximations are better than those obtained 
from performing the eigenvalue analysis with only 4 blade elements. 
Furthermore reference 7 suggests an iteration method to improve the 
accuracy of any specific eigenvalue, that was found to rapidly converge 
to the exact value in a numerical example. From Table 5 there seems to 
be no doubt that for the rot or- airframe coupling problem modalization of 
the polynomials with the single hinged blade eigenvalues and subsequent 
truncation of the partial fraction expansion is a viable method to reduce 
the order of the total system characteristic equation, if this is required. 
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Concluding Remarks 

In the oonventlonal stabiUtyanalysis for coupled rotor-alrframe 
systems structural modes are used where all parts of the system 
oscillate in phase. In order to approximate the actual modes, where 
due to aerodynamic effects substantial phase differences exist between 
the oscillations of the varioiis parts of the system, several structural 
modes must be superimposed, and it is difficult to judge how many of the 
structural modes are needed to properly represent a true aeroelastlc mode. 
In the finite element stability analysis explored here for the case of 
coupled blade flap-bending and airframe oscillations the aerodynamic 
effects are directly incorporated into the blade element transfer 
matrices. The resulting natural modes are true aeroelastic modes and 
the problem of structural mode superposition is eliminated. Computational 
limitations occur due to numerical errors from repeated multiplications 
of high order polynomials. For 10 blade elements and performing the 
computations with 16 digits the numeiical errors present no problem. 

The analysed system is then of 44th order. Though the method should be 
tried out on more sophisticated descriptions of both the rotor and the 
airframe to include blade chordwise bending and blade torsion and a more 
complex airframe structure, it appears from the experience gained so far, 
that the principal limitation is in the largest eigenvalue considered. 

It can be anticipat^id that with more eigenvalues but approximately 
retaining the largest on^ no numerical difficulties will be encountered, 
except that the order of the characteristic equation may become too high 
for accurate root extraction. In this case the method of modalizing and 
truncating the polynomials, briefly studied herein, could be applied to 
reduce the order of the total system characteristic equation. 
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While the main purpose of this investigation was to shed somt 
light on the tractability of the finite element st utility analysis for 
coupled rotor- airframe systems, interesting results were obtained in 
the numerical examples, [or the selected ratio of blade flapping 
moment of inertia over rotor support inertia of 5, and for the selected 
hingeless blade first natural flap-ben din frequency vntiiout damping 
and while rotating of 1.06 there are 4 str '.’gly coupled aeroelactic 
modes, 2 modes that are for small and for large support stiffness 
mainly support modes, and 2 further modes that are in these regions 
mainly progressing first blade- flap bending and regressing second 
blade-flap bending modes. Two of these modes can become unstaijle or 
near unstable in the intermediate range of support stiffness. The 
progressing second blade flap bending and the higher blade benaing 
modes show little tendency of coupling with the rotor support and 
retain high damping. The extension of the finite element stal>ilitv 
analysis to forward flight conditions should pose no great diffienJi ie.. 
as long as the rotor advance ratio is moderate and constant cocfficier.t 
in the multiblade equations can be considered adequate. 
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